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ABSTRACT 

We explore the linear stability of astrophysical discs exhibiting vertical shear, which 
arises when there is a radial variation in the temperature or entropy. Such discs are 
subject to a “vertical-shear instability”, which recent nonlinear simulations have shown 
to drive hydrodynamic activity in the MRI-stable regions of protoplanet ary discs. We 
first revisit locally isothermal discs using the quasi-global reduced model derived by 
Nelson et al. (2013). This analysis is then extended to global axisymmetric perturba¬ 
tions in a cylindrical domain. We also derive and study a reduced model describing 
discs with power law radial entropy profiles (“locally polytropic discs”), which are 
somewhat more realistic in that they possess physical (as opposed to numerical) sur¬ 
faces. In all cases the fastest growing modes have very short wavelengths and are 
localised at the disc surfaces (if present), where the vertical shear is maximal. An 
additional class of modestly growing vertically global body modes is excited, corre¬ 
sponding to destabilised classical inertial waves (“r-modes”). We discuss the properties 
of both types of modes, and stress that those that grow fastest occur on the shortest 
available length scales (determined either by the numerical grid or the physical viscous 
length). This ill-posedness makes simulations of the instability difficult to interpret. 
We end with some brief speculation on the nonlinear saturation and resulting angular 
momentum transport. 

Key words: accretion, accretion discs - planetary systems - hydrodynamics - waves 
- instabilities 


1 INTRODUCTION 

Accretion through magnetorotational turbulence is only vi¬ 
able in sufficiently ionised regions of protoplanetary discs, 
namely at their inner and outer radii ( [Balbus Sz Hawl^ 
[1998 Armitag^|2011 ). Between 1-20 AU (the “dead zone”) 
non-ideal effects extinguish the MRI, and instead accre¬ 
tion may occur via magneto-centrifugally launched winds 
(e.g. Lesur et al.||2014] Bai||2014 ). However, identifying ad¬ 
ditional hydrodynamic mechanisms for driving turbulence 
is essential, due to its potential impact on the dynamics of 
solids, and therefore for planet formation. 

Though pure Keplerian shear flow is difficult to desta¬ 
bilise, several mechanisms have been proposed: subcritical 


paloizou 2010), convective instability (|Ruden et al. 1988 


baroclinic instability (Petersen et al. |2007 Lesur & Pa- 


1964 


Lesur Ogilvie|2010 ) and gravitational instability (Toomre 


Lin & Pringle 1987), to name but a few. Another mech¬ 


anism that has recently received attention is the “vertical- 
shear instability” (hereafter VSI) which, as its name sug¬ 
gests, attacks rotating systems that exhibit vertical shear 
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( Urpin Brandenburgp998 Urpin|[2003 ). Fundamentally, 
the VSI is a form of centrifugal instability and is a close 
cousin of the Goldreich-Schubert-Fricke instability, origi¬ 


nally applied to stellar interiors (Goldreich & Schubert 1967 


Fricke||1968[). But obse rvations of protostellar discs ( An¬ 
drews V Williams|2005 ) and theoretical models of passively 
heated discs ( Ghiang &; Goldreich||199'7 ) suggest that they 
too should display destabilising vertical shear, generated 
from radial variations in temperature or entropy. Recent nu¬ 
merical simulations indicate that the nonlinear evolution of 
the VSI can produce hydrodynamic turbulence and modest 
levels of angular momentum transport ( Nelson et al.|[20T3 


Stoll V Kley|2014 ). It thus could be a potentially key player 


in the dynamics of protoplanet ary disc dead zones. 

The VSI was originally studied with a local (Boussi- 
nesq) approach by Urpin V Brandenburg (1998) and by 
Urpin| ( [2003 ). More recently Nelson et al. ( 2013| described 
it with a quasi-global model that captures the full vertical 
structure of growing anelastic modes in radial geostrophic 
balance (assuming the background is locally isothermal). Be¬ 
ing global in the vertical, but local in the radial, this model is 
akin to the commonly used vertically stratified shearing box. 
However, several properties of the linear VSI require further 
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explanation, especially with respect to its global manifes¬ 
tation in more realistic disc models. This is a particularly 
important issue when trying to connect the linear theory to 
global simulations, and in interpreting their nonlinear out¬ 
come. Our paper is devoted to exploring this aspect of the 
problem. 

We perform linear stability analyses of astrophysical 
discs exhibiting global variations in temperature and en¬ 
tropy, and as a consequence vertical shear. We employ lo¬ 
cally isothermal and polytropic models in both quasi-global 
and fully global 2D geometries, which revise and extend pre¬ 
vious work. 

In agreement with Nelson et al. (2013), we find that 
the VSI excites two classes of modes. The first class cor¬ 
responds to classical free inertial waves (r-modes) that are 


present in any astrophysical disc (Lubow & Pringle 


1993 


Korycansky & Pringle 1995 Kato 2001) but which have 


been destabilised by the vertical shear. These, referred to 
as “body modes”, grow at modest rates and typically ex¬ 
hibit longer wavelengths (though the radial wavelength of 
the waves is still short). 

The second class corresponds to modes localised to the 
vertical surfaces of the disc where the vertical shear is max¬ 
imal. These grow much faster and have very short wave¬ 
lengths, making them difficult to resolve numerically. In fact, 
unless viscosity is included, the fastest growing modes pos¬ 
sess arbitrarily small wavelengths, making their simulation 
problematic. Note that, though they have been termed “sur¬ 
face modes”, these are different to the classical surface grav¬ 
ity waves that appear in poly tropic disc models, as they lie 
in a different frequency range; they are hence a form of lo¬ 
calised low-frequency inertial wave. Strict isothermal mod¬ 
els do not possess a physical vertical surface and hence do 
not support these surface modes. Polytropic disc models do, 
however, as should any realistic disc model that possesses a 
transition between an optically thick interior and an opti¬ 
cally thin “corona”. 

We begin by explaining why a radial variation in en¬ 
tropy or temperature generally leads to vertical shear in § [^ 
There we also explain why such discs are likely to be unsta¬ 
ble. After defining our basic disc models in § [^ we analyse 
the resulting VSI in the locally isothermal disc in §|^and[^ 
and the locally poly tropic disc in § [^Finally, we will discuss 
the implications of our results in where we also specu¬ 
late on the nonlinear evolution of the VSI and its efficiency 
at transporting angular momentum. 


2 VERTICAL-SHEAR INSTABILITY 

Discs with radial variations in temperature or entropy neces¬ 
sarily possess vertical shear. To see that this must be, con¬ 
sider the “thermal wind equation” (the azimuthal compo¬ 
nent of the vorticity equation for the axisymmetric basic 
state of the disc): 

= -e^ • (Vp X VP)/p" (1) 

= duTd^s - d.TdRS. ( 2 ) 

Here we have adopted cylindrical polar coordinates centred 
on the central object (R, 0, z) and p, P, S and T are the 
basic state density, pressure, specific entropy and tempera¬ 
ture profiles, respectively. Eq. tells us that a radial varia¬ 


tion in the background temperature or entropy generates a 
departure from cylindrical rotation through the baroclinic 
terms on the right hand side. Thus the angular velocity 
Q, = D(P, z), and consequently the disc exhibits a weak 
vertical shear. For illustration, we show the angular velocity 
and vertical shear for a disc with a radial variation in tem¬ 
perature in Fig. and the vertical shear for a disc with a 
radial variation in entropy in Fig. (both disc models and 
the notation adopted are defined in § |^ . 


2.1 Physical picture 


Vertical shear provides a source of free energy that can drive 
hydrodynamic instabilities. How might modes access this 
free energy? Consider a ring of fluid at a given location 
(A) within the disc with coordinates {Ra,za), and hence 
specific angular momentum Ha = R\Q{Ra, za)- If we ver¬ 
tically displace this ring to a new position (B) with coor¬ 
dinates {Ra,za + Az), then its specific angular momen¬ 
tum will be conserved as long as viscosity can be neglected 
(i.e. we assume that jAzj is much larger than the viscous 
length). But if the angular momentum of fluid at the new 
location hs is smaller (larger) than Ha, then the ring will 
be pushed outwards (inwards) by the centrifugal acceleration 
{h\ — h%)/R\, leading to a dynamical instability. Given that 
h% ^ h\ Azdzh^, instability occurs whenever dzh? < 0 
(or indeed > 0), i.e. if there is any vertical shear. Basically, 
this interchange of rings of fluid reduces the total energy of 
the configuration, leading to an instability that transports 
angular momentum in order to eliminate the vertical shear 
This is a modified form of Rayleigh’s argument for centrifu¬ 
gal instability. Though accretion discs are stable according 
to the classical Rayleigh criterion, any vertical shear permits 
its circumvention and hence the onset of instability. 

This physical argument works for neutrally stratified 
discs, but must be altered when a stable vertical entropy 
stratification is present, as exhibited by protoplanet ary disc 
dead zones. So we next introduce buoyancy, and viscous and 
thermal diffusion. Buoyancy forces impede exchanges of the 
type described, and thus inhibit any adiabatic (dynamical) 
instability (cf. the Solberg-Hpiland criterion). Instability is 
nevertheless possible if the buoyancy forces are eliminated, 
such as by sufficiently fast cooling or thermal diffusion. For 
this to work, displacements |Az| must then be much shorter 
than the thermal diffusion scale. The resulting instability is 
hence double-diffusive in character, possessing lengthscales 
lying in the range bounded from below by the viscous length 
and above by the thermal diffusion length. Originally identi¬ 


fied in the 1960’s and applied to stellar interiors (Goldreich 
fc Schubert 1 1967 Fricke|19^ ), only much later was it recog¬ 


nised that such an instability could emerge in astrophysical 
discs ( Urpin V Brandenburg||l998| )p] 


^ Our illustrative perturbation is vertical for simplicity; any dis¬ 
placement lying within the angle between the rotation axis and 
a surface of constant angular momentum will do (as explained in 
K nobloch &: Spruit ] 1982' for ex ample). 

^ I Urpin ^ Brandenburg ( |l998| coined the term “vertical shear 
instability”. However, there is a good case for the retention of 
the name “Goldreich-Schubert-Fricke” (GSF) instability, even if 
“VSI” has the merit of clearly advertising the underlying physics. 
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(a) a (b) d^(m) (c) p 


Figure 1. Basic state for the locally isothermal disc with q = — l,p = —1.5 and cq = 0.05. The left panel shows a contour plot of 
on the (R, z) plane. The middle panel is a similar contour plot, but this shows the magnitude of the vertical shear dz{R^), which has a 
maximum at \z\ ~ 1 (whereas the scaleheight at the inner radial boundary is 0.05). The right panel shows the density p. 


2.2 Estimates from a local model 


According to a local Boussinesq analysis, the growth rate of 
the VSI is 




(3) 


where e — H/R\s the disc aspect ratio and q is the exponent 
in the power law for temperature (or entropy), so that T oc 


R^ ( Urpin Brandenburg|p^^ Urpin|[2003) [Nelson et al. 


2013). For \q\ ~ 1, the VSI will hence grow, and presumably 
saturate, relatively quickly on a timescale not far from the 
dynamical time for thicker discs. 

The VSI afflicts intermediate lengthscales i in the range 


< i < 4 


(4) 


where the viscous and thermal diffusion lengths are defined 
through 

e, = {,yla)K 4 = (x/iV.)h (5) 


Here ly is the kinematic viscosity, x is the thermal diffusiv- 
ity, and Nz > 0 is the vertical buoyancy frequency. On the 
other hand, lengthscales above are stabilised by buoyancy 
forces (as long as > 0); only when these are subdued by 
sufficiently rapid thermal diffusion is instability possible. On 
lengthscales smaller than , viscosity neutralises the excess 
angular momentum of a fluid element too quickly to allow 
the mode to grow. The fastest growing modes have vertical 
(kz) and radial wavenumbers (kn) that satisfy 


kz 

kn 




eq, 


( 6 ) 


SO the radial wavelength of the mode (in) is typically much 
shorter than the vertical wavelength (iz) ( Urpin V Bran- 
denburg|1998 Urpin|2003 ). Hence, being inertial waves, the 
group velocity points almost vertically, in accordance with 
the direction of transport outlined in the previous subsec¬ 
tion. Note that modes grow at or near the fastest rate on 
short lengthscales all the way to the viscous cut-off 

To give some sense of the numbers, the growth rate for 
the VSI is typically a ~ 0.3 yr“^, if e ~ 0.05 at 1 AU and 
\q\ ~ 1. The microscopic kinematic viscosity at the mid¬ 
plane of a protoplanetary disc at 1 AU may be estimated as 
iy rsj 2.5 X lO^cm^s”^, yielding ~ 10^ km, much shorter 


than the local disc thickness H (of the order of 10^ km). 
The thermal diffusivity, on the other hand, is significantly 
larger, x ~ 5 x lO^^cm^s”^ (taking values typical for H 2 in 
a Minimum Mass Solar Nebula). Note that x varies signif¬ 
icantly with height in the disc (e.g. [Bell Lin 1994), but 
we ignore this complication when making simple estimates. 
It is unclear what value Nz should take in a dead zone 
in a protoplanetary disc, but it is probably much smaller 
than Q. Therefore a (very crude) lower bound on the ther¬ 
mal diffusion length is ^ 10^ km. In all likelihood, how¬ 
ever, H, not least because thermal diffusion becomes 

stronger as we approach the photosphere. 

These estimates suggest the following ordering: H ~ 
^ with the VSI occurring in the wide gulf separating 
the viscous and thermal diffusion lengths. We might expect 
the fastest growing modes to localise in regions of greatest 
vertical shear; in any realistic disc model, the magnitude 
of the vertical shear increases with distance from the mid¬ 
plane, and takes its maximum value at the disc surfac^ So 
we might expect the fastest growing modes to occur on very 
short lengthscales (just above h, ~ 100 km) located near the 
surface. 


2.3 Nonlinear evolution 

We expect that the VSI will work to eradicate the destabil¬ 
ising conditions from which it arose (the vertical shear) and 
ultimately return the system to a marginally stable, cylindri¬ 
cal rotation prohle. Of course, resisting the VSI will be the 
driver of the vertical shear itself: the radiation held of the 
protostar. In the struggle between these two opponents, the 
system will probably reach a quasi-steady ‘balance’ in which 
the vertical shear is diminished (but not entirely removed) 
and some degree of hydrodynamical activity remains. The 
properties of this state will be determined by the relative 
efficiency of the VSI in wiping away the shear versus the 
shear’s thermal driving via the protostar. Presumably, if the 

^ This is not the case for the locally isothermal thin disc model 
discussed in §3.1.1| which has no surface and for which the mag¬ 
nitude of the vertical shear takes it maximum value at vertical 
infinity, where the density is negligible. 
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Figure 2. Basic state for the locally polytropic disc with qs = 
— l,p = —1.5, ao = 0.05 and 7 = 1.4 (therefore H{R = 1 ) ps 0.13). 
The top panel shows a contour plot of the vertical shear dz (R^) 
on the (i^, 2 ;)-plane. For given R^ the magnitude of the vertical 
shear has a maximum near the disc surface. The bottom panel 
shows the density p. Regions outside the disc surface where p = 0 
are coloured white in both panels. (The angular velocity is not 
shown since the departure of constant surfaces from the vertical 
is very small). 


VSI is inefficient then significant vertical shear persists and, 
as a consequence, significant turbulent motions|jThe locally 
isothermal simulations of Nelson et al. (2013) best describe 
this scenario, because the destabilising gradients are fixed 
and cannot be modified by the VSI. On the other hand, if 
the thermal driving is weak, the VSI should eliminate the 
vertical shear, and subsequently its motions will settle down 
to a much lower level, because the system is near (if not at) 
marginal stability. This situation can be observed in simu¬ 
lations where the unstable equilibria possess long relaxation 
times ( Nelson et al.||2013 ) or in cases where the disc is not 
thermally driven ( Stoll fc Kley||2014 ). 

However it is generated, the properties of the final quasi¬ 
steady state are of key importance to dead-zones, interven¬ 
ing potentially in the dynamics of solids and even in angu¬ 


lar momentum transport. It certainly should be the focus 
of future numerical simulations. Note that this state will 
evolve slowly, on the long timescale of the disc’s and proto¬ 
star’s evolutionary track, and also presumably on the shorter 
timescale of the protostar’s emission variability. 

The VSI in protoplanet ary discs is challenging to sim¬ 
ulate adequately because the fastest growing modes have 
very short lengthscales, many orders of magnitude shorter 
than both the disc thickness and (typically) the numerical 
grid scale. Global simulations are then not only ill-posed but 
exhibit a nonlinear saturation whose characteristic length- 
scales are forced to be longer than is realistic. It may be 
that global simulations of the VSI greatly overestimate the 
amount of power in the largest lengthscales. Indeed [Stoll 
Kley (2014) find no convergence in radial angular momen¬ 


tum transport (a) as resolution is increased, a result that 
emphasises the prominence of the smallest scales. 

We can crudely estimate an upper limit for the result¬ 
ing radial angular momentum transport if we assume that 
~ \q\^^z and if nonlinear saturation occurs when the ve¬ 
locity amplitude is of order Irg. Consequently, the turbulent 
viscosity scales as rt ^ ^rCt ~ \q\^e^ To make further 
progress we must estimate the size of iz- If the dominant 
scales at saturation are the largest ones, then £z H, which 
gives an upper bound on the a parameter: 


< 




( 7 ) 


For a protoplanet ary disc with e 0.05, we get a < 10 
which is bro adly consistent with the results of curr ent global 


simulations (Nelson et al. 


2013 


Stoll fcKley|2014| |p1lf, how¬ 


ever, we assume instead that the dominant scale at satura¬ 
tion is much shorter (so that £z H), then this is a gross 
overestimate. 


3 HYDRODYNAMIC EQUATIONS AND 
BASIC STATE PROFILES 

In the previous section, we explained why discs with global 
radial variations in temperature or entrop^j^are likely to be 
unstable to the VSI. We now list the equations and describe 
the simplified disc models that will be used to analyse the 
VSI. We begin with the equations of compressible hydrody¬ 
namics for an inviscid adiabatic fluid: 

(at + M-V)m=--VP-V4>, (8) 

p 

(dt U • V) p = -pV • u, (9) 

(ft + 'a-V)>S = 0, (10) 

where u is the velocity. The ideal equation of state is P = 
IZpT (where IZ is the gas constant divided by the mean 
molecular weight). 

If we adopt cylindrical polar coordinates (P, 0, z), the 


^ In this case, the fastest growing VSI-driven modes probably 
saturate through secondary shear instabilities - as is the case for 
fingering convection, another kind of double-diffusive instability 
( [Brown et al.|[^13| ). The resulting turbulence in the low den¬ 
sity regions near the disc surface may best be studied using local 
simulations, since they can access more realistic small scales. 


^ Note that many of the a calculations in jStoll &: Kleyj ( |2014| 
are two dimensional and questions may be asked of angular mo¬ 
mentum transport in this case (see arguments in Balbus|2000 ). 

® We do not set out to analyse the stability of local structures 
in the thermal properties of the disc, such as edges or pressure 
bumps etc. 
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gravitational potential due to the central object is approxi¬ 
mately that of a point-mass 




GM 


V-R2 + 22 ’ 


GM 


1 - 




= $0 + $22 , 


( 11 ) 

( 12 ) 


where in some cases we expand for a thin-disc (|«| <C R', 
second line), and we define $o = —GM/R and $2 = 
GM/{2R^). 


3.1 Basic state 


The axisymmetric basic state of the differentially rotating 
disc has u = RQ(R, z)e(i), and satisfies the equations of ra¬ 
dial and vertical force balance: 

-Ra^en = --VP + V$. (13) 

P 


Taking the curl then provides the thermal wind equation 
(Eq.0. For simplicity, and to allow some analytical reduc¬ 
tion and exploration, we restrict ourselves to studying two 
simple models (as in Nelson et al.|2013 ): the locally isother¬ 
mal disc with a radial power law in temperature, and the 
locally poly tropic disc with a power law entropy function. 
In both cases, we want to consider a disc with mid-plane 
density 


Pm ('^) 



(14) 


where po is the mid-plane density at a radius Rq. It turns 
out that the value of p is unimportant for the VSL 


3.1.1 Locally isothermal disc with a radial power law in 
temperature 


The first model that we will discuss is a locally isothermal 
disc with P = (?s{R)p (where Cs is the isothermal sound 
speed), in which the temperature depends only on cylindri¬ 
cal radius as the power law 

T{R) = n(^^y, (15) 

and similarly c^s{R) — Cq{R/RqY, where Cq = IZTq is the 
square of the isothermal sound speed at a radius The 
corresponding density profile that satisfies Eq. |13| is 

p{R, z) = pm{R) exp [MR) - «)]) • (16) 

Note that the disc has no surface and formally extends to 
infinity. The angular velocity that satisfies Eq. is 

2 \ i 

niR,z) = MR)(i + q + ip + q)§^- ’ (i^) 


where QoiR) = V^M/R^ and H{R) = c,{R)/QoiR) oc 
(R/Rop'^^''^'^ is the local disc scaleheight. The angular ve- 
locity therefore depends on z whenever q Y 0. We set out 
to consider p and q so that these discs are stable according 
to the Solberg-Hpiland criteria governing adiabatic axisym¬ 
metric perturbations (e.g. Tassoul||1978 ). 

Eor illustration, in Eig. we plot the angular velocity 


Q, vertical shear dz{RLl) and density p on the {R, z) plane 
for a disc with cq = 0.05,p = —1.5 and q = —1. This shows 
that the angular velocity depends on z, and that the verti¬ 
cal shear increases monotonically with z until its reaches a 
maximum far away from the mid-plane, where the density 
is negligible. 


3.1.2 Locally polytropic disc with a radial power law in 
entropy 


The second model is a locally polytropic model in which P — 
Ks{R)p^ (where 7 is the adiabatic index) with an entropy 
function 

Ks{R) = Pp--^ = Ko (18) 

with qs and Kq constants. Hence Ks is only a function of 
cylindrical radius, so that S = S{R) oc \nKs{R). This disc 
model is neutrally stratified in the vertical direction. The 
corresponding density profile which satisfies Eq. p^is 

piR,z) = p^{R)y + Yy^[MR)-m,z)]ym 

where m = 1/(7 — 1) is the polytropic index. Eor the last 
line, the potential has been expanded for a thin disc, and we 
have defined the local disc thickness 


Ho{R) = 


2(1 + m)Ks{R) 


( 21 ) 


where Qo was defined in § |3.1.1| This disc model possesses 
a surface at which p = 0 when z = Ho{R). The adia¬ 
batic sound speed is a{R, z) = which 

becomes am{R) at the mid-plane (taking the value ao at 
R = Ro), and we define M(R,z) = RQq{R)/ a{R, z) and 
Mm{R) = RLlo{R)/am{R) to be the Mach number and mid¬ 
plane Mach number of the flow. The angular velocity that 
satisfies Eq. [^is 


niR,z)-noiR){i + qs + ^ + :^-^^== 


■ ( 22 ) 


As earlier, the angular velocity depends on z whenever 
7 ^ 0. Since this model is neutrally stratified in the vertical 
direction, it can become unstable to adiabatic disturbances 
for certain choices of p and qs , whenever one of the Solberg- 
Hpiland criteria are violated. However, we do not set out to 
study such instabilities in this work, and we instead focus 
on discs that would be stable to adiabatic perturbations. 

Eor illustration, in Eig. we plot the vertical shear 
dz{RLl) and density p on the {R^z) plane for a disc with 
ao = 0.05,7 = = —1.5 and qs = —1. This shows that 

the angular velocity depends on z, and that the magnitude 
of the vertical shear increases monotonically with |z| until 
its reaches a maximum just below the disc surface. 

In the next three sections we work through the stability 
of these two equilibria. Because the full global analysis is 
challenging, we first treat their quasi-local approximations 
in a reduced model first outlined in Nelson et al. (2013). 
This helps tease out the most important features. Once this 
is done we perform the full global analysis on the locally 
isothermal model only. 
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LINEAR STABILITY OF THE LOCALLY 
ISOTHERMAL DISC: REDUCED MODEL 


In this section we revisit the reduced model of INelson et al.l 
(2013), which describes the dynamics of the VSI in locally 


isothermal discs with a radial power law in temperature. 
We begin by outlining the derivation of this model (a similar 
derivation is presented in detail in Appendixj^for the locally 
polytropic disc) and go on to analyse its most important 
properties. 

We define e — HjR and consider a thin disc in which 
e <C 1. We are interested in slow modes with frequencies 
and growth rates that are each O(eQ), with vertical scales tz 
that are comparable with the thickness of the disc [tz ~ eR) 
and radial scales Ir that are much smaller {tn ~ e^R), cf. 
Eq.H If we also consider vertical velocities that are mildly 
subsonic or transonic (by a factor e) and radial velocities 
that are very subsonic (by a factor e^), then we require 0(e) 
density perturbations, for consistency. On such small radial 
scales, the curvature of the disc can be neglected and the 
geometry is locally Cartesian, similar to the classical shear¬ 
ing box ( Goldreich Lynden-Bell|1965 Umurhan & Regev 


|2004| . In this limit, we are looking at low frequency (iner¬ 
tial) dynamics that are anelastic and in approximate radial 
geostrophic balance. 

The linearised reduced equations for the rescaled veloc¬ 
ity perturbation (i^, w) and the fractio nal density pertur¬ 
bation Ii — p/~p are ( Nelson et al.|2013 ) 

(23) 

(24) 


0 

drV 

drW 

0 


2v - dxU, 
u qzw 
~2 

-a.n. 


= dx(pu) + dz(pw). 


(25) 

(26) 


Here r is rescaled time, and x is a local radial variable. The 

^2 

background density is p = e“^, after appropriate normali¬ 
sation. The crucial term for the appearance of the VSI is the 
last one in Eq. |24[ which arises from thermal wind balance 
if there is radial variation in temperature (see Eq. [^. 

We seek solutions of this system of the form 


n = 


= Re 


(27) 


and so on for other variables, where we subsequently drop 
the hats for clarity. This allows us to reduce the system to 
a single ODE that can be written most simply by defining 
new (complex) coordinates ( = z^/]A^Lq, as 

d^n ^dH ^ 

^-C^ + AR-O, 

where we have set 
cu k 


A = 


1 + ikq' 


probabilist’s version), which has solutions 

n = ai HeA(C)+a 2 

The first function is the Hermite function and the second 
function is a confluent hypergeometric function of the first 
kind, with ai, a 2 arbitrary constants. Solutions that are 
polynomially bounded as |C| ^ oo (so that \p'\ 0 as 


Id ^ oo) require A = n G N and a 2 = 0 (see Okazaki 
|et al.||1987| |Kato||2001 ). The regular solutions are therefore 


nocHen(C), (31) 

with Hcn a Hermite polynomial of order n. The correspond¬ 


ing uj can be obtained from the dispersion relation Eq. (29). 
These solutions describe the vertical structure of modes in 
the low frequency limit. Since they are (complex) polynomi¬ 
als in d they describe global modes which are not localised 
in the vicinity of any particular location C 7^ 0- 


4.1 Non-vertically shearing case, g = 0 

Before treating the VSI we examine the case when g = 0 
in order to make contact with the existing literature on 
wave modes in vertically stratified isothermal discs (namely 
Lubow V Pringle|1993 ). This then can assure us of the valid¬ 
ity of the reduced model and the regime of its applicability. 

Since A is quantised, we obtain a discrete set of frequen¬ 
cies. When g = 0, these are real and 

j_Vn 


(32) 


which represents a pair of low frequency inertial wave^trav- 
elling in opposite directions. The full isothermal disc allows 
axisymmetric waves with frequencies that satisfy the follow¬ 
ing dispersion relation ( [Lubow V Pringle|1993 ) 

(—+ 1) — (?sk^ uA — 0. (33) 

There are two branches of solutions, which represent ei¬ 
ther high frequency acoustic waves or low frequency inertial 
waves. Note that there are no surface modes in the isother¬ 
mal disc, since the model lacks a surface. The inertial waves 
have ~ n/{lA + n), therefore Eq. 32 is a good approxi¬ 
mation for modes with k ^ as expected from the as¬ 
sumption IrAz ~ e. When this is violated, physically incor¬ 
rect solutions appear with frequencies larger than one. The 
restriction on the size of the vertical quantum number n 
means that the model can never capture purely local modes 
(for which both n and k are large). 


4.2 Vertically shearing case, ^ / 0 

When g 7 ^ 0 the disc is unstable to the VSI. The frequencies 
are now 


(28) 

o; = ±A^(l + ifcg)L 

(34) 


which are complex. The growth rate a 
computed from the positive root of 

= Im[ct;] may be 

(29) 


(35) 

(the 

The unstable modes in this case are the classic inertial waves 
of the vertically stratified disc, destabilised by the vertical 

(30) 

shear. They are all global “body modes” 
isation near to any C 7^ 0- 

that have no local- 


^ Note that these are sometimes referred to as r-modes in the 
literature (and sometimes inappropriately as g-modes even when 
buoyancy forces are absent). 
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Vertical-shear instability 7 


In the limit \kq\ ^ 1, 

0 J = ±^ I 1 + iy+0(fc q ) 


(36) 


which possesses the same frequency as the classical modes, 
cementing the identification. They grow, however, at the 
rate a ^ ^/n\q\/2. In the opposite limit \kq\ ^ 1, 




(37) 


therefore the growth rate is cr = ^n|g|/( 2 /c), which is no 
longer linear in the shear. 

Note that the growth rate increases with polynomial 
order as ^/n, which may seem somewhat surprising but can 
be understood by considering the vertical shear profile. An 
isothermal thin-disc lacks a surface, yet oc \qz/2\, 

which increases without bound as |z| ^ oo. Since modes 
with large n extend over a greater vertical range, they can 
better tap into the energy associated with the larger shear 
at larger z. Hence it follows that the growth rate should in¬ 
crease with n (the square root comes from the fact that n 
has units of inverse length squared). It should be remem¬ 
bered, however, that these expressions apply only when the 
reduced model is valid, which requires ^/n k. 

Finally, note that the unstable modes here have a wave 
character, whereas those issuing from the local Boussinesq 


calculation grow monotonically (Urpin & Brandenburg 1998 


Urpin 2003). Their scalings with wavenumber also differ. 


In the isothermal limit with k/kz ^ 1 and k ^ local 
Boussinesq modes grow at a rate a ~ ^/k^qJk, where kz is 
the vertical wavenumber. The system, being scale free, only 
depends on the ratio of wavenumbers, kz/k. In contrast, 
the quasi-global model separates out the vertical and radial 
scales, forcing the former to be near H. Consequently, it 
will always struggle to reproduce the local limit. Only if we 
force modes to localise at a fixed z (such as at a boundary) 
and then undertake a plane-wave analysis, can the reduced 
model produce the local Boussinesq results ([Nelson et al.| 

Ml. 


4.3 Vertically shearing case with imposed 

boundaries: appearance of surface modes 

To connect with recent numerical simulations of a locally 
isothermal disc ( Nelson et al.||2013 Stoll k, Kley||2014" ), we 
adopt an artificial boundary at a finite height in the ver¬ 
tical direction. In realistic discs there is surely some form 
of transition at the photosphere, but whether a numerical 
boundary adequately mimics this feature is unclear. We em¬ 
phasise that within the confines of a strict isothermal model, 
a numerical boundary is an artificial addition, yet it has im¬ 
portant effects. 

We consider boundary conditions such that the verti¬ 
cal momentum p^w = 0 at |z| = i7, where H is free pa¬ 
rameter (other choices would lead to similar behaviour for 
the VSI). To illustrate the unstable modes in this case, we 
solve Eqs. [23H^ numerically using a Chebyshev collocation 
method on V +1 points of a Gauss-Lobatto grid. This results 
in a matrix eigenvalue problem that can be solved using a 
QZ method ( Boyd| 200 T] Golub van Loan|1996 ). Numeri¬ 
cal convergence for the eigenvalues is verified by varying the 
resolution and comparing eigenvalues. 



Figure 4. Growth rate of the unstable modes versus their (real) 
frequencies in an isothermal disc with artificial boundaries at 
\z\ = if, illustrating the lack of convergence as H is varied for 
cases with A; = 10 and q = —1 (converged results were obtained 
with N = 300). Note that the maximum growth rate in this case 
is approximately \qH\/2. We have also indicated the analytical 
prediction of Eq.j^with the solid red line to guide the eye (how¬ 
ever, it should be remembered that this represents a set of discrete 
modes and not a continuum). The branch that extends approxi¬ 
mately vertically represents surface modes. 


In the presence of an artificial boundary, a new class 
of “surface modes” appears that localise near to the bound¬ 
aries. These are in addition to the “body modes” obtained 
previously. We show the vertical momenta for several exam¬ 
ples of each type of mode in Fig. for k = 10 and q = —1: 
the top panel shows a selection of typical surface modes, 
and the bottom panel shows typical body modes. The mode 
plotted in the top right panel is one of the slowest growing 
surface modes, and its character is intermediate between the 
two classes of mode. The appearance of these two classes of 


modes was also found by Nelson et al. (2013) and McNally 


& Pessah (2014). 


Mathematically, the emergence of surface modes can 
be understood by examining Eq. |30| In the presence of a 
boundary at a finite height, A is complex and ai and a 2 are 
both nonzero, in general. New modes, localised near to |z| = 
i/, appear due to the need to match boundary conditions 
at this location, rather than at infinity (this can be verified 
by plotting HeA(C) for A G C). Thus the new surface modes 
rely entirely on the boundary for their existence. 

In Eig. we show the dependence of the spectrum on 
disc height H. Here the numerical growth rates of the un¬ 
stable modes are plotted in the complex frequency plane for 
three different heights, H = 5, 6 , 7. The remaining parame¬ 
ters are k = 10 and q = —1. We also plot the if = oo case 
Eq. |34| which is given by the solid red line in the figure. 
Roughly, surface modes correspond to the more “vertical” 
segment of the spectrum, and the inertial waves to the more 
“horizontal” segment. 

Obviously, there is no convergence with increasing if. 
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Figure 3. Illustration of the real and imaginary parts of the vertical momenta for a representative selection of the two types of modes 
in an isothermal disc with imposed numerical boundaries. Here k = 10, q = —1, H = 5 and N = 300. The figure labels indicate the 
complex frequencies of the modes. Surface modes are shown in the top row (the fastest growing mode is in the top left panel), and body 
modes are shown in the bottom row. Similar surface modes localised at the bottom boundary are also obtained. The lowest frequency 
body mode is plotted in the bottom right panel and is an n = 1 (n = 0 for w) inertial wave, i.e. the fundamental “corrugation mode”, 
which is well described by an n = 1 {n = 0 for w) Hermite polynomial in the absence of vertical boundaries, and its complex frequency 
is that predicted by Eq. |34| 


especially for the surface modes. As the height of the domain 
is increased, the maximum growth rate increases in direct 
proportion with H. This is what we would expect, since 
<7max oc max|^ 2 (i^Q)| = \q\H/2 (even larger growth rates 
than shown in Fig. [^are obtained for modes with k ^ 10, 
which asymptotically attain this value). The presence of a 
boundary also strongly influences the inertial waves, only 
the lowest n of which are well described by Eq. As H 
is increased, however, modes with increasingly larger n con¬ 
verge to the if = oc analytical prediction. 

The appearance and lack of convergence of the surface 
modes is a special pathology of the isothermal model with 
imposed boundaries. It makes the interpretation of these 
surface modes, and their role in any ensuing turbulence, es¬ 
pecially problematic. Are they merely numerical artefacts? 
In the next section we argue that they are, in fact, more than 
that and that it is the vertically isothermal model itself that 
is the problem. 


5 LINEAR STABILITY OF THE LOCALLY 
POLYTROPIC DISC: REDUCED MODEL 

In this section we analyse a reduced model of the VSI in discs 
with a radial power law in entropy. The vertical structure 
of this locally polytropic model is a good approximation for 
an optically thick disc, and is more realistic than the locally 
isothermal model because it possesses upper and lower sur¬ 


faces. In addition, there is a well-defined maximum vertical 
shear rate, which occurs at the disc surface. Therefore this 
model is better defined mathematically and physically. 

A formal derivation of the reduced model is relegated to 
Appendix The same scalings as in § |^ are adopted, cor¬ 
responding to anelastic radially geostrophic phenomena. In 
addition, thermal diffusion is likely to dominate the thermo¬ 
dynamics of such slow and short-scale modes. Consequently, 
we may neglect the entropy perturbation entirely, as well as 
the stabilising influence of stratification. 

The resulting linearised equations for such perturba- 


tions are 



0 = 

2v - dxh, 

(38) 

drV = 

u qsZ 

~2 

(39) 

II 

-dzh, 

(40) 

0 = 

„ , „ 2mzw 

OxU + dzW 2 ’ 

(41) 


where the pseudo-enthalpy perturbation is h = P'/{I — 
and P' is the (scaled) pressure perturbation. Rapid thermal 
diffusion means that the perturbations evolve isothermally, 
hence the equations are similar in form to Eqs 123126] 

We seek solutions of the form 

u = Re , (42) 

and so on for other variables. The resulting system of equa- 
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Figure 5. Comparison of the dispersion relation computed nu¬ 
merically from the locally polytropic reduced model (Eqs. |38| - 
1 ^ with gfg = 0 (dashed lines) together with the numerically 
computed predictions for inertial modes in a polytropic disc ( |Ko-| 
[rycansky &: Pringle] 1995t (solid lines) with n = 1, 2, 3, 4, 5 vertical 
nodes adopting 7 = 1.4. The reduced model is only valid when 

k > 0(yn). 


tions is equivalent to the single ODE 


(Vh 

d«2 




2m ikqs 

1 — z‘^ 7 


27 27 _ n 

—— T ^ k h — 0. 
dz 


(43) 


Unfortunately, this cannot be solved in closed form analyti¬ 
cally. However, we can attack the problem using a matched 
WKBJ approach, or numerically with the same method as 
§ |4.3| Chebyshev collocation on A^ + 1 points, followed by the 
eigensolution of the resulting matrix equation using the QZ 
method. No explicit boundary conditions are imposed, but 
implicit regularity is assumed at |z| = 1 (this gives identical 
results to imposing a free-surface condition explicitly). 

Finally, we note that the resulting system can be writ¬ 
ten in a scale-free manner by transforming our variables to 
hatted quantities, as follows: 


V = kv, w = kw, uj = uj/k, Qs = 2q^lk. (44) 


A similar rescaling also applies to the isothermal model. Un¬ 
der this rescaling the governing equation depends solely on 
the two parameters q and m. Each, discrete, growth rate 
then has the following behaviour: 


= j^LOn{kqs, j), 


(45) 


where the function Un is determined numerically. 


5.1 Non-vertically shearing case, = 0 

We first demonstrate that the reduced model correctly cap¬ 
tures the inertial waves in a polytropic disc in the limit of 
large k, in order to gain confidence in the reduced model. 

The full polytropic disc model without vertical shear 
{qs = 0 ) possesses three different classes of neutrally stable 


axisymmetric modes: high-frequency acoustic modes, low- 
frequency inertial modes, and surface gravity modes of in¬ 
termediate frequency ( Korycansky Pringle|p^^ Ogilvie 


1998). The reduced model is designed to capture only the 


low frequency inertial modes. 

In Fig. we compare the frequencies of five inertial 
modes with vertical modenumbers (vertical nodes) n — 
1,2, 3,4, 5 and = 0 for several k as predicted by the re¬ 
duced model and the full polytropic model of |Korycansky 


Pringle (1995) (solving their Eqs. 4-8). The latter are rep¬ 


resented by solid and the former by dashed lines. Unlike 
the isothermal disc, these frequencies must generally be ob¬ 
tained numerically. There is general agreement for /c > 20, 
which illustrates that the reduced model correctly captures 
the frequencies of these inertial modes for sufficiently large 
k. However, only modes with k ^ 0(v5r) are correctly cap¬ 
tured. 

In the limit of large k, a WKBJ analysis shows that the 
frequencies take a simple form. 


Ct; = 7 r( 2 n + m)/(4/c), (46) 

where n is an integer. Though the formula is most accurate 
in the limit of large n, it does well across the whole range 
of frequencies. Details of this calculation can be found in 
Appendix B. 


5.2 Vertically shearing case, qs ^ 0 

When qs ^ 0, the disc is unstable to the VSI. As an illustra¬ 
tive example, we plot in Fig. the complex frequencies of 
the unstable modes computed from the eigenvalue problem 
Eqs. [38}{4T] when k = 100, qs = —1 and 7 = 1.4 for two 
different vertical resolutions N = 200 and N = 400. The 
unstable modes (for this k) are shown to be well resolved 
using N = 200. 

Fig.| 6 ] also shows clearly that a locally polytropic disc 
can be unstable to two different types of modes: (a) modestly 
growing inertial waves (“body modes”; these occur on longer 
radial scales, as we will show in Fig.|^, and (b) rapidly grow¬ 
ing short-wavelength surface modes (these only occur when 
^ 30, as we will also show in Fig.[^. The growth rates 
of the two classes of modes are labelled in Fig. Examples 
for each type of mode are plotted in Fig. Hence we find 
that the VSI in the locally poly tropic disc, with well-defined 
physical surfaces, is qualitatively similar to the VSI in the 
isothermal disc with an imposed artifical surface at a finite 
height, as discussed in § |4.3| 

In Fig. we plot the complex frequencies of the un¬ 
stable modes when — 1 and 7 = 1.4 for several values 

of k. To accurately capture the fastest growing modes for 
larger k requires increasing the vertical resolution {N = 400 
was required to obtain convergence when k = 300). Fig. 
shows that the growth rate of the fastest growing mode in¬ 
creases, and its frequency decreases, as we consider smaller 
horizontal lengthscales (larger k). For a given |^s|, the num¬ 
ber of rapidly growing surface modes increases with k. All 
unstable modes are body modes when k = 10 , but as we 
increase k surface modes appear, and by /c = 300 there 
is a large population of them. The maximum growth rate 
is on track to approach the maximum vertical shear rate 
|^ 2 (i^fl)| ~ \qs\H/{2'y) ^ 0.357 as /c ^ oc, which we expect 
to provide an upper bound on the growth rate of the VSI. 
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Figure 7. The real and imaginary parts of the vertical momenta for several representative examples of both types of modes in a polytropic 
disc with k = 100, = — 1, 7 = 1.4 using N = 200, where the figure labels indicate the complex frequencies of the modes. In the top 

row we show several examples of surface modes, with the fastest growing modes plotted in the two top left panels. In the bottom row 
we show several examples of body modes. The lowest frequency mode is plotted in the bottom right panel, and is an n = 1 (n = 0 for 
w) inertial wave, i.e. the fundamental “corrugation mode”. The other bodes modes are higher order body modes. 
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Figure 6. Growth rates of the unstable modes versus their (real) 
frequencies in a polytropic disc with k = 100 , qs = — 1 , and 
7 = 1.4 computed using N = 200 and N = 400 points. This 
shows that the modes are well captured using N = 200, and the 
distribution is similar to the isothermal case in Fig.^ 


Figure 8. Growth rates of the unstable modes versus their (real) 
frequencies for several values of in a polytropic disc with qs = 
— 1 and 7 = 1.4, computed using up to AT = 400 points (which was 
found to be sufficient for all k). This illustrates the dependence 
of the unstable modes on k. 


In Fig. we have exploited the rescaling of Eq. 
to illustrate the general dependence of the scaled growth 
rate cf/IqsI of the VSI as a function of k\qs\ for all unstable 


modes. We do not consider modes with real frequencies that 
are larger than unity - their omission results in a region 
where the VSI is absent in the bottom left of the figure. 
Pairs of surface modes appear for sufficiently large k\qs\. 
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Figure 9. Normalised growth rates of all unstable modes as a 
function of the scaled radial wavenumber k\qs \ in a polytropic disc 
with 7 = 1.4, computed using N = 200. The absence of unstable 
modes in the lower left of the figure results from our elimination 
of modes with frequencies larger than 1 . Surface modes appear 
in pairs for k\qs\ > 30, with the fastest growing arising from 
the n = 2 “breathing mode” and the n = 3 inertial mode. The 
n = 1 fundamental “corrugation mode” lies along the line that 
ends up at the bottom right. The scaled growth rate of the fastest 
growing surface mode asymptotically approaches the maximum 
vertical shear rate (0.357) as k\qs \ oo. 


The first, and fastest growing, pair of surface modes appear 
when k\qs\ > 30, and corresponds to the n = 2 “breathing 
mode” and the n — 3 inertial mode changing their charac¬ 
ter. The n — 1 fundamental “corrugation mode” is the first 
to become unstable when k\qs\ >3, but its growth becomes 
weaker for larger k\qs\. This is represented by the curve that 
crosses all others and ends up at the bottom right of the fig¬ 
ure. As k\qs\ oo, there are more unstable modes, and their 
scaled growth rates approach the scaled maximum vertical 
shear rate 1 /( 27 ) ~ 0.357, as expected. 

The fastest growing body mode is the n = 1 corru¬ 
gation mode, and has a maximum growth rate of approxi¬ 
mately 0 . 1 |gs|( 1 . 4 / 7 ) when \kqs\ < 10. The fastest growing 
mode is a surface mode, with a maximum growth rate that 
approaches the maximum vertical shear rate |^ 2 (i^Q)| ~ 

I 1 /( 27 ) on the smallest scales. (We have confirmed the de¬ 
pendence of the fastest growing surface mode on 7 , though 
we omit this for brevity - this arises because the basic disc 
structure varies with 7 .) Note that the full polytropic disc 
possesses a class of surface modes that are restored by grav¬ 
ity (Korycansky & Pringle 1995 Ogilvie 1998). However, 
these are unrelated to the class of low frequency surface 
modes that we have discussed in this section, and occur in 
a different frequency range. 

Why do surface modes appear when |A;^s| > 30? The 
vertical shear rate takes its maximum magnitude at the disc 
surface, so it makes sense that if localised modes were to 
appear, they should do so just below the disc surface. But 
when \kqs\ < 30, it is not possible for a mode (with an iner¬ 


tial wave character) to become sufficiently localised close to 
the disc surface (i.e. iz cannot be made sufficiently small). 
The criterion for the appearance of surface modes is there¬ 
fore tied to the shortest vertical lengthscale permitted. 

Finally, in Appendix B we show that the body modes’ 
growth rate can be obtained analytically in the limit of small 
Qs: 


mqs 


7 r 7 ( 2 n + m] 


log [| 7 r( 2 n + m)] , 


(47) 


where n is a (large) integer. The corresponding wave fre¬ 
quency is given by Eq. It should be stressed that the 
growth rate a being at subdominant order is only a rough 
estimate, because errors arising from the WKBJ method it¬ 
self enter at the same order. Note that the growth rate is 
linear in the shear, but the larger n, the smaller cr, in con¬ 
trast to the isothermal case, but in agreement with Figs. 
and For general qs an estimate for the surface modes’ 
maximum growth rate is cr ~ ^ 5 /( 27 ), while its wave fre¬ 
quency scales as Ink/k. Being leading order, this is a more 
robust estimate, which is in approximate agreement with 
our results. 

In summary, the VSI in the locally polytropic model is 
very similar in character to the VSI in the locally isothermal 


model with artificially imposed vertical boundaries (§ 4.3). 
The presence of a surface, be it physically justified or a nu¬ 
merical convenience, provides a special location upon which 
modes can affix themselves and localise. The polytropic disc 
however yields solutions that are clearer to interpret because 
the location of the surface is specified (not an adjustable 
parameter). Consequently, there exists a well-defined maxi¬ 
mum growth rate for the instability, set by the vertical shear 
rate at the disc surface. 


6 GLOBAL CALCULATIONS IN THE 
LOCALLY ISOTHERMAL DISC 

We have so far analysed the VSI using reduced models ol 
locally isothermal and polytropic discs. Here we present the 
first two-dimensional stability calculations ol the VSI in a 
locally isothermal disc. One motivation lor doing so is to 
verily the validity ol the model analysed in § |^ another 
is to reproduce the instability in a setup that more directly 
matches that ol recent global simulations ( Nelson et al.|20f^ 
Stoll Kley|[20T4 ). An issue with global stability calcula¬ 
tions is that they are computationally demanding (primarily 
regarding their memory usage), so we are limited to study¬ 
ing the VSI with relatively low resolutions. This is primarily 
a problem lor capturing the surlace modes, since they occur 
on very short lengthscales. However, the lowest order body 
modes have less vertical structure and are better captured 
in these calculations. 

We consider axisymmetric perturbations {ur, U(j),Uz, p') 
to the global disc model ol Section 3.1.1, assuming their 
time-dependence to be oc . The equations governing 
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their linear evolution are then 


- lUJUR = 2RQu^ - dR{clp) + ^dR{clp) + fR, (48) 


-iuu^ = -UR — dR{R^Q) - Uzdz{R^) + , 
R 

-iujuz = -^dz{Llp) + ^dz{clp) + fz, 


(49) 

(50) 


-lujp = -urOrp - Uzdzp - ^dR{RuR) - pdzUz,(51) 

rC 

where P' = cl{R)p', and the background state quantities 
Cs(R),^{R, z) and p{R,z) are defined in § 3.1.1 To regu¬ 
larise the solutions in some calculations we include a Navier- 
Stokes shear viscosity with constant kinematic viscosity u, 
with the extra terms 


f = 


V^u + -V(V • It) 


+ vp'V'^Uo 


+iyVp ■ |^V« + (Vm)^ - -(V • m)I 
+rVp ■ [vC/o + (Vt/o)^] , 


(52) 


where C/o = RQ^e^. We adopt a 2D cylindrical domain with 
R G and z G [—zo,^o], where zq is some multiple 

of i/o = e (this differs from the spherical wedge considered 
by [Nelson et al. (2013) and [Stoll Kley (2014), but this 
difference is probably unimportant). Our units of length and 
time are chosen such that Ro = 1 and Q(i?o) = 1. If / = 0, 
boundary conditions are enforced such that 


puR = 0 at R = Ro & R = Ri, 

puz = 0 at z = —zo & z — zq^ 


(53) 

(54) 


otherwise we supplement these with stress-free conditions 
(no tangential stresses) on each boundary. This system is 
solved (after writing the equations in terms of momenta 
rather than velocities) numerically using a Chebyshev collo¬ 
cation method in both R and z, with Nr-\-1 and A/’^-fl points 
on a Gauss-Lobatto grid, respectively. The resulting gener¬ 
alised eigenvalue problem involves matrices of size L x L, 
where L — A{Nr + l)(iV 2 H- 1), and is solved using one of two 
methods: a QZ algorithm to obtain an approximation of the 
full spectrum, or an Arnold! iterative method to obtain an 
approximation to a desired part of the spectrum. The QZ 
method is computationally expensive when L > 20,000, so 
our results are then supplemented by the Arnold! method, 
once we determine which modes to focus on. In what follows 
we set p = 0, since preliminary investigation suggested this 
parameter to be unimportant. 

The code has been tested in several ways. First, without 
viscosity and with a small radial domain (e.g. Ri = 1.001), 
we have confirmed that our code accurately reproduces the 
inertial and acoustic modes of the vertically unbounded 
isothermal disc when g = 0 ( Lubow fc Pringl^|1993 ) - ex¬ 
cept for modes with large n, for which the confining effect 
of the vertical boundaries modifies the solutions. Using a 
similarly small radial domain, we have verified that viscous 
damping produces the expected decay rate for a mode with 
a given short radial wavelength {vk^ when 1). In all cal¬ 
culations, the discretisation inevitably produces many junk 
modes, that involve oscillations on the grid-scale. Where 
possible, we eliminate these modes by comparing eigenvalues 
for several different resolutions, and by inspecting the spa¬ 


tial structure of the eigenfunctions, discarding those that 
oscillate on the grid-scale. 


6.1 Inviscid calculations 

We first illustrate the properties of the VSI in the absence 
of viscosity. As in the reduced model in § ^ we obtain two 
classes of modes: modestly growing body modes, and rapidly 
growing surface modes localised near the numerical bound¬ 
aries m. z. As in Section 4.3, these modes would vanish if the 
boundaries could be taken to infinity. However, we do not 
discard these modes entirely as they still probably represent 
physical solutions in some sense (see § . 


6.1.1 Body modes 


As in § |4.2| we obtain a set of body modes, essentially classi¬ 
cal inertial waves that grow in the presence of vertical shear. 
We plot the vertical momenta for several examples with ei¬ 
ther even or odd symmetry in Fig. in a domain with 
= 1.5 and zq = 5e. These modes have vertical structures 
with n = 1,2,3 and 4 (n = 0,1,2,3 for pw), which corre¬ 
spond with those obtained in the isothermal model without 
boundaries in § 4.2 (and the bottom panels in Fig.[^ for var¬ 
ious radial structures. The vertical shear is strongest at the 
inner radial boundary of the domain, whereas these modes 
are concentrated near the outer boundary. We might ex¬ 
pect this to be the case because axisymmetric inertial waves 
are localised at and within their turning surfaces, defined 
by K = cj (where r ^ Q is the epicyclic frequency). The 
best resolved modes in our calculations have low frequencies, 
and thus their turning surfaces lie near or beyond the outer 
boundary Ri. Consequently such modes prefer the largest 
radii possible in the computational domain. Those plotted 
here were chosen to illustrate that the global model exhibits 
body modes with the same vertical structure as the reduced 
model in § |4.2[ The body modes are modestly growing, with 
a growth rate no larger than a third of the maximum vertical 
shear rate. These are the unstable modes with the longest 
wavelengths, and have a radial scale that is shorter than 
their vertical scale by a factor 0(e). 


6.1.2 Surface modes 

Just as in § |4.3| we obtain a class of rapidly growing short- 
wavelength surface modes, which come in pairs with oppo¬ 
sitely signed frequencies. For illustration, we show an exam¬ 
ple of a typical surface mode with and without viscosity in 

Fig.[T| 

Without viscosity, the fastest growing mode always oc¬ 
curs on the smallest available lengthscale. To demonstrate 
this, in the top panel of Fig. we plot the growth rate and 
(real) frequencies of the fastest growing mode versus resolu¬ 
tion iV in a disc with e = 0.05 and g = — 1, adopting a do¬ 
main with Ro = 1.0025 for several different resolutions N = 
Nr = g [20, 30, 40, 50, 60]. Such a small radial domain is 
chosen in order to better capture the rapidly growing sur¬ 
face modes. We also plot the vertical momenta for the fastest 
growing mode for several different resolutions in Fig. |13| As 
we increase the resolution, the fastest growing mode moves 
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(c) OJ = 0.2019 + 0.0034i (d) oj = 0.2209 + 0.0035i 


Figure 10. Illustration of the vertical momenta on the (R, 2 :)-plane for several body modes in the locally isothermal disc without viscosity 
with e = 0.05, q = —1^ zq = and Ri = 1.5, obtained using Nji = Nz = 40. These modes have vertical structures for Uz that are well 
described by Hermite polynomials with n = 0,1,2 and 3. Many modes with similar vertical structures (but different radial structures) 
to each of these are also unstable. Also over-plotted are three contours of constant density (solid black lines) with p = 10 “^, 10 “^, 10 “^ 
for reference. The amplitude in each panel is arbitrary. 


to lower frequency and exhibits increasingly shorter length- 
scales. In addition, its growth rate increases as we increase 
the resolution, gradually tending towards the maximum ver¬ 
tical shear rate (which is |g|i7o^o/(2i?o) ~ 0.125). This is 
because these modes become increasingly localised in the 
vicinity of the vertical boundary, where the vertical shear is 
maximal. Note that Fig.[^does not indicate convergence as 
N is increased, because it is a different mode that is plotted 
(that is most unstable) for each iV, and this mode always 
tracks the grid-scale. This is demonstrated further in Fig.p^ 
Since these modes always occur on the grid-scale, they are 
necessarily the most poorly resolved modes. 

Another pathology of the isothermal model, that we 
first showed in § |4.3[ is that the maximum growth rate de¬ 
pends on the vertical domain size. To further demonstrate 
this, we illustrate the fastest growing mode for calculations 
with zo = He, with H G [5, 6, 7, 8] on Fig. in a radial 
domain with Rq = 1.0025. We plot results using both a 
fixed resolution of Nr = 20 and Nz = 60 (blue circles and 
line), and for fixed Nr = 20 and NzjH — 12 (black crosses 


and line). The growth rate continues to increase as we in¬ 
crease H. This is simply explained by the dependence of the 
vertical shear rate, which increases without bound as 0{z). 
However, the finite vertical resolution does not fully cap¬ 
ture the fastest growing mode for each H (even with fixed 
Nz/H), therefore the dependence on H is slightly weaker 
than the linear extrapolation (red dashed line). 


The maximum growth rate for the surface modes com¬ 
pares reasonably well with the numerical simulations of |Nel-| 


son et al. (2013) for similar parameters (they obtain a growth 


rate of 0.125 whereas our maximum growth rate is approx¬ 
imately 0.1 for N = 60, for example - the difference is due 
to the inevitably smaller resolution in our case). However, 
growth rates are unconverged with respect to both (a) the 
numerical grid resolution and (b) the vertical domain size. 
The fastest growing mode occurs on the smallest available 
lengthscale, which is always the grid-scale. In addition, the 
fastest growing mode occurs at the boundaries in z. Since 
the boundary is artificial imposed and this choice sets the 
maximum growth rate, our results are strongly dependent 
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(a) o) = 0.0257 + 0.0757i (Nr = = 20) (b) oj = 0.0174 + 0.0893i (Nr = 


30) (c) oj = 0.0126 + 0.0965i (Nr = = 40) 


Figure 13. Illustration of the vertical momenta on the {R, z) plane for the fastest growing surface mode without viscosity for several 
different resolutions, in a locally isothermal disc with e = 0.05, q = —1, zq = 5e and Ri = 1.0025. These have increasing resolution from 
left to right panels, as indicated in the figure labels. This shows that the fastest growing surface modes occur on the smallest available 
lengthscales, hence this problem is ill-posed because modes on the numerical grid scale will always be the fastest growing. Note the 
eompressed radial seale, which was chosen to capture such rapidly growing surface modes - these modes are actually inclined by a small 
angle to the vertical, as we expect. The amplitude in each panel is arbitrary. 



Figure 11. Growth rate and frequency of the fastest growing 
surface mode as a function of numerical resolution N = Nji = Nz 
in an isothermal disc without viscosity, with e = 0.05, q = —1, 
Zq = 5e and Ri = 1.0025. These results were computed using 
an Arnoldi method. The growth rate increases with resolution N, 
and the frequency correspondingly decreases, indicating a lack of 
convergence. The red-dashed line is the maximum vertical shear 
rate. 


on an arbitrary parameter H, which far from ideal. These 
problems will also afflict nonlinear simulations of the VSI in 
locally isothermal models. 



(a) Nr = Nz = 50,1/ = 10-8,a; = -0.0419 + 
0.068H 



(b) Nr = 30, Nz = 120, u = 0, LB = 0.0092 + 
0.0607i 


6.2 Viscous calculations 

Although the dependence on the vertical boundary cannot 
be removed, convergence with respect to resolution can be 
addressed via the inclusion of viscosity. This is necessary for 
the problem to be well-posed. As we will demonstrate here. 


Figure 12. Illustration of the vertical momentum on the {R, z) 
plane for a typical surface mode with and without viscosity, in an 
isothermal disc with e = 0.05, q = —1, zq = 5e and Ri = 1.01. 
The amplitude is arbitrary. 
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(a) Ro = 1.5, Nii = = 60 


(b) Ro = 1.05, Nu = 30, Nz = 150 


Figure 15. Growth rate of unstable modes versus their (real) frequencies in a locally isothermal disc with e = 0.05, zq = 5e, Ri = 1.5 
and 1.05 for several viscosities (resolution is listed in the figure labels). The dot-dashed line in the left panel demarcates the low-frequency 
VSl modes from the viscously overstable modes. The VSl appears only when u < 10“®. As the viscosity is decreased, there are more 
unstable modes, which occur on shorter lengthscales, and their maximum growth rates increase. (We caution that some of the eigenvalues 
in the right panel for u = 10 “® in particular do not appear to be converged with the adopted resolution - however, we do not expect 
any of our conclusions to be adversely affected by this.) 



Figure 14. Growth rate of the fastest growing surface modes 
versus vertical domain size in a locally isothermal disc without 
viscosity for several values oi zq = He, with e = 0.05, q = —1 and 
Ri = 1.0025. Results were computed using an Arnoldi method 
with fixed resolution Nji = 20 , Nz = 60 (blue circles and line) 
and with Nji = 20 and a fixed value of Nz/H = 12 (black crosses 
and line). The growth rate of the fastest growing mode increases 
with H, indicating the lack of convergence as we increase the 
vertical size of our domain. 


even when we include viscosity, the dominant modes prefer 
the smallest available lengthscales. 

In Fig. |15| we plot the unstable modes on the complex 
frequency plane for converged modes in a domain with zq = 
5e in various radial domains with several values of u. Note 
that these figures are busier than say Fig. because they 
include modes of all resolvable vertical and radial quantum 


numbers; in Fig.[^ the radial wavenumber k is restricted to 
take only one value. 


The left panel shows unstable modes when Ro = 1.5. 
Two classes of instability are obtained in this case: the VSl 
with uj < 0.4 and the viscous overstability with cu > 0.4, 
which have been separated visually by the blue dot-dashed 


vertical line. The viscous overstability (Kato 


1978 


Kley 


et al. 1993| [Latter V Ogilvie 2006) preferentially excites 


modes with little vertical structure and long radial wave¬ 
lengths k~^ > 0{Hq), with frequencies comparable with the 
local rotational frequency and growth rates 0{iyk^). This 
can be distinguished from the VSl, which occurs only when 
q ^ 0, preferentially excites waves with short radial wave¬ 
lengths, Ct; ^ 1, and has growth rates bounded above by 
the maximum vertical shear rate. When g = 0, only the vis¬ 
cous overstability persists - this is modified by vertical shear 
when q = —1, but is found to remain in the same region of 
the complex plane. Given that our focus here is on the VSl, 
we will not discuss these modes any further. 


The VSl is absent if the viscosity is too large {u > 10“®) 
since it preferentially excites (in this case body) modes with 
short radial wavelengths which are damped by a weak vis¬ 
cosity. We note that a similar order of magnitude for u is 


required to obtain the VSl in numerical simulations (Nelson 
et al.||2013 ). The fastest growing mode when iv = 10“ is a 


body mode that has n = 1 (i.e., the fundamental “corruga¬ 
tion mode”) with a radial wavelength of approximately 0.05 
(similar to the top left panel of Fig. 10). 


The right panel in Fig. plots the unstable modes in 
a smaller radial domain of Ro = 1.05, for several viscosities. 
The smaller domain is chosen so as to permit a greater radial 
resolution. This demonstrates that as we decrease the vis¬ 
cosity, there are (a) more unstable modes and (b) the fastest 
growing mode has a larger growth rate, in tune with our ex¬ 
pectations from the reduced model. The increase in the num- 
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ber of unstable modes is because the smaller v the shorter 
the viscous cutoff; consequently, there is a wider range of 
lengthscales that are potentially unstable (cf. Fig.|^. When 
V — there is only one unstable mode in this case, which 
we omit for clarity. 

We also observe that the fastest growing modes occur on 
the shortest available lengthscales, since the mode frequency 
decreases with decreasing v. However, we note that there 
are many unstable modes, with the fastest growth rate be¬ 
ing somewhat smaller than the maximum magnitude of the 
vertical shear (0.125) as a result of viscous damping. Surface 
modes appear, and become dominant when ia < 10“^. 


6.3 Summary 


Our results in this section are in accord with the reduced 
model introduced by Nelson et al. (2013) and revisited in §[^ 
but with imposed vertical boundaries. We obtain the same 
classes of modes and they have similar growth rates in both 
models (and compared with numerical simulations) in a fi¬ 
nite vertical domain. We have highlighted that the VSI pref¬ 
erentially excites ultra short-scale disturbances which occur 
on the smallest available lengthscales, be they numerical or 
viscous. Adoption of explicit viscosity is therefore required 
to obtain results that are converged with respect to resolu¬ 
tion. 


True isothermal discs do not exhibit surface modes, as 
we have explained in §[^ Only the presence of a rigid bound¬ 
ary permits those modes to exist, but such a model is poorly 
defined given the freedom regarding our choice of zq. Given 
the notable lack of convergence as the vertical domain size 
is varied, this indicates that the vertically isothermal model 
is not well suited for studying the linear properties of the 
VSI. 


7 DISCUSSION AND CONCLUSIONS 


We have analysed the linear stability of astrophysical discs 
with vertical shear, as a result of their radial variations in 
temperature or entropy. Such variations are expected to be 
present in real discs, as indicated by both observations and 
theory (e.g. Andrews V Williams|2005' Chiang V Goldreich 
1997), and generally lead to vertical shear, thereby rendering 


the disc unstable to a hydrodynamic instability. Recent non¬ 
linear simulations of the resulting vertical-shear instability 
have highlighted its potential to drive hydrodynamic activity 


in MRI-stable regions of protoplanetary discs (Arlt & Urpin 


2004 Nelson et al. 2013 Stoll V Kley 2014). The aim of 


this work was to better understand the nature of the linear 
instability in two simple disc models: the locally isothermal 
disc with a radial power law in temperature (building on 
previous work by Nelson et al.|[20T3 ) as well as the locally 
polytropic disc with a radial power law in entropy. 

In both models, there are two classes of unstable 
modes: modestly growing (vertically) global body modes, 
and rapidly growing ultra short-scale surface mode^ The 


® The surface modes are somewhat analogous to the modes that 
appear in fingering convection, another type of double-diffusive 
instability (e.g. Brown et al.|2013 ). 


latter only appear in discs with a vertical surface and, 
though this is not the case in strictly isothermal models, re¬ 
alistic discs should exhibit a density feature/transition upon 
which such modes can affix themselves. Ironically, artificially 
imposed boundaries present qualitatively correct behaviour 
even if this behaviour is a numerical artefact! The value of 
such models beyond the qualitative is unclear, however. 

A separate issue is that surface modes preferentially oc¬ 
cur on the smallest available length scales. This necessitates 
the inclusion of viscosity to obtain convergence in any nu¬ 
merical simulation, otherwise the results will inevitably de¬ 
pend on numerical resolution. It may be that the VSI sat¬ 
urates on length scales much smaller than can be reached 
by nonlinear global simulations, which would consequently 
pump power to artificially larger scales and hence misrepre¬ 
sent the ensuing turbulent state. This however will be very 
difficult to test numerically. 

We have restricted our study to the locally isothermal 
and locally polytropic models, since studying the linear VSI 
in a more realistic model would require the two-dimensional 
disc structure to be computed numerically after accounting 
for the various sources of heating and cooling, which are un¬ 
certain. In any case, it is unlikely that there is much to be 
gained from doing this, owing to the similarity in its prop¬ 
erties in both the locally isothermal disc (with an artificial 
boundary) and the locally poly tropic disc. 

As shown by recent work, the nonlinear evolution of the 
VSI leads to wave activity or turbulence, which transports 
angular momentum vertically in order to eliminate the ver¬ 
tical shear ( Arlt V Urpin|[2004 Nelson et al.||2013) Stoll V 


Kley 2014). It may also transport angular momentum ra¬ 


dially to enable the disc to accrete at modest levels. The 
amount of such hydrodynamical activity in reality depends 
on the battle between the external heating of the disc, and 
the efficiency of the VSI in eliminating the ensuing verti¬ 
cal shear. The disc must be externally heated sufficiently 
strongly or the VSI will eventually win out, leading to a 


very low level or no turbulence (as observed in Stoll & Kley 


2014). That this process can occur on timescales that are 


not very long may preclude the use of local computational 
models to determine the transport properties of the VSI. 
This is unfortunate, because global simulations of the VSI 
that are able to capture the fastest growing modes (even in 
the presence of viscosity) are computationally very challeng¬ 
ing. Nevertheless, it would be worthwhile to determine the 
nonlinear outcome of the VSI in more realistic disc mod¬ 
els (continuing from Stoll & Kley 2014) to determine its 
longevity and transport properties. 

An interesting byproduct of the VSI’s saturation is its 
radial transportation of angular momentum. If the VSI is 
sufficiently active in protoplanetary dead zones, it may pro¬ 
vide the stresses necessary for these regions to accrete. The 
saturation of the VSI could be controlled by the smallest- 
scale surface modes, which may be ineffectual in this regard 
because of their small scales. On the other hand, longer 
wavelength body modes carry greater quantities of angu¬ 
lar momentum, but only if the saturated state endows them 
with sufficient power. If this is the case a crude upper limit 
on the a associated with this transport is of order ~ 10“^ 
(which is roughly consistent with previous simulations). On 
the other hand, Stoll V Kley (2014) show that the measured 
a decreases with resolution, and is moreover unconverged. 
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Vertical-shear instability 17 


a result that suggests that it is in fact the smallest scales 
that are controlling the transport, not the largest, and that 
in real discs the value of a could be negligible. Obviously 
much more work is required to further test these ideas. 
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APPENDIX A: DERIVATION OF A REDUCED 
MODEL FOR THE VERTICAL-SHEAR 
INSTABILITY IN A LOCALLY POLYTROPIC 
DISC 


We look for linear axisymmetric 
{uR,Uct>,Uz,P',p',S') to Eqs. H 


10 


perturbations 
to which we also 


include thermal diffusion in the entropy equation of the 
form • [xVT], where x = x{p,T,k) = is the 

thermal conductivity. We will end up neglecting thermal 
diffusion when constructing the basic state, but will include 
it when considering the small-scale perturbations. The mo¬ 
tivation is that we are considering slow perturbations with 
very short radial lengthscales, on which thermal diffusion 
is very rapid, so the perturbations evolve approximately 
isothermally. 

The system of linearised equations is 


dtUR = 

2mu,p - -OrP' + ^OrP, 

p 

(Al) 

dtUfj) — 

-UR^dRiR^Q) - u,dz{RQ), 

JrC 

(A2) 

dtUz = 

--dzP' P ^dzP, 

p 

(A3) 

dtp = 

—URdRp — Uzdzp — ^dR(RuR) — pdzUz, 
rC 

(A4) 

dtS' = 

-urOrS - u,dzS + ^ V • (xVT') . 

(A5) 


We adopt a non-dimensionalisation such that our fidu¬ 
cial radius has Rq = 1 and Qo = 1, and take e = Hq/Rq = 
Hq as our small parameter. To obtain a reduced model, we 
consider slow dynamics on vertical scales comparable with 
the disc thickness and radial scales that are much smaller. 
We also assume that the resulting velocities are comparable 
with the sound speed, with the exception of the radial ve¬ 
locity, which is assumed to be much slower. In particular, 
we define a slow timescale r — et such that dt — edr, along 
with new radial coordinate x and vertical coordinate z, such 
that 


R - 1 = e^x, Or = e ‘^dx, (A6) 

Z = ez, dz — e~^dz- (A7) 


This allows us to neglect curvature and to consider the ge¬ 
ometry as locally Cartesian, similar to the classical shear¬ 


ing box ( Goldreich fc Lynden-Bell|1965 Umurhan & Regev 


2004). We also define scaled velocity components (u,u,ie) 


such that 

Ur = U(f) = eu, Uz — ew, (^8) 

and appropriately scale the perturbed pressure, density and 
entropy for the approximation to be consistent as follows: 

P' = VP, p = tp, s' = eS. (A9) 

The background pressure is also rescaled as P = e^Po- We 
also note that Or acting on background quantities (e.g. p) is 
0(1) and dz acting on background quantities is 0(e“^). Un¬ 
der the above scaling assumptions (and neglecting thermal 
diffusion), the basic state in § 3.1.2| can be written 

2 1 m 

1 - , (AlO) 


Po 


P = Pa 


Hi 


2 1 m+1 


1 - 


R^} — RqQ^q 


Hi 
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where the disc thickness is 


Ho = J2{l + m) — . 

Po 


(All) 

(A12) 

(A13) 


Finally, we take x = xV, where /3 is an ordering pa¬ 
rameter. We require /3 > 2 so that thermal diffusion can be 
neglected when constructing the basic state, since 


- pT{u -VS) = V • (xVT) 

= [xdiT + d.xdzT] 

+e<^ [xdiT + dRxdRT] , 


(A14) 


where the left hand side is 0(1). For the perturbations, we 
have 


V • (xVT') = 


e-^+^XdiT + e-^^^XdiT 
^~^^^dRxd^f + t-^+^d,xdzf. 


-1+P,^2r, 


+e“ 


(A15) 


The dominant term here is clearly the first one, and this 
dominates over all other terms in Eq. |A5| as long as /3 < 3, 
since the leading order term on the left hand side of Eq. |A5| 
is wdzS, which is 0(1). If /3 < 3, the influence of a stabilis¬ 
ing entropy gradient on the perturbations is eliminated by 
rapid thermal diffusion. We will therefore consider {3 G (2, 3). 
While this choice may seem contrived, what this means 
physically is that we are neglecting thermal diffusion for 
the basic state, but we are considering it to dominate the 
thermal evolution of the perturbations, which evolve isother¬ 
mally. 

Applying these scaling assumptions to Eq. |A1}{A^ leads 
to the reduced model 


0 

II 

to 

1 

(A16) 

drV 

_ U QsZ 

(A17) 

drW 

= -dzh 

(A18) 

0 

= p [dxU + dzw] + wdzP 

(A19) 


at leading order, where the psuedo-enthalpy perturbation 
is h = P/P- The perturbations are in radial geostrophic 
balance and are anelastic. 
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APPENDIX B: WKBJ ANALYSIS OF THE 
LOCALLY POLYTROPIC REDUCED MODEL 

In this section we obtain approximate analytic solutions to 
Eq. 1^ in the limit of large k. To ease the calculation, we 
first rescale u and Qs , introducing the quantities u = kuj and 
q = kqs/{2^). The limit of large k then becomes the limit of 
large to and q. Next Eq. [^is transformed into Schrodinger 
form 


dz^ 


^2 1 — (m — l)z^ 


(1-z2)2 

I .. 1 — (1 + 2m)z‘^ ^2 2 

+iq -^ 


if = 0, (Bl) 


for the new dependent variable H = ^‘^h{z). 


Note that asymptotically Eq. |B1| is the harmonic oscillator 
equation for all z except for small regions near the bound¬ 
aries at z = ±1. 

Near the upper boundary z — Eq. |B1| may be ap¬ 
proximated by 

d^H 


dz^ 


+ 


'.2 , .2 , m{2-m) 
UJ + q + 


imq 

1-z 


H = 0, (B2) 


4{l-zA 

which admits a solution in terms of generalised Laguerre 
functions and the Tricomi hypergeometric function. Only 
the former, however is regular at the boundary. The appro¬ 
priate solution is hence 

if w (1 - - z) 

r m—1 


(B3) 


where ^ is the generalised Laguerre function, zj = 


\J—CS^ — and 


m f iq 


(B4) 


In the limit of large lu, the solution takes the following 
asymptotic form far from the boundary, 

H ~ r(?Ti - fi) exp[0(2:)] + r(fi) exp[-0(2:)], (B5) 

where T is the gamma function, and the phase function is 

ITTfl 


= 117(1 — z) -\- IBtAL ln[ 2 ii 7 (l — z) 
2w 


+ (B6) 


This solution, in fact, should hold throughout the rest of the 
domain and so we impose boundary conditions (evenness or 
oddness) at the mid-plane. This yields the eigenvalue equa¬ 
tion 

r(M) 


^7 + log(2ra7) - I log 


where n is an integer. 


r(m — jj) 


+ |i7r(n + /i) = 0, 
(B7) 


Equation (B7) is transcendental in w. However, in the 
non-vertically shearing case it can be solved easily and we 
obtain the very simple dispersion relation 


(jO — — {2n + m ). 
4/c 


(B8) 


recalling that n is an integer, and m is the polytropic index. 

Eor general g, Eq. |B7| must be solved via a root-finding 
algorithm, usually an easier task than tackling the ODE 
itself. However, when ^ ~ 1 (meaning kqs <C 1), we find 
that the wave frequencies of the unstable body modes are 
given by Eq. |B8| to leading order, while the growth rate is 
mqs 


7T'y(2n + m 


■ log [|7r(2n + m)] . 


(B9) 


The growth rate is linear in the shear, but the larger n, the 
smaller a, in contrast to the isothermal case. This therefore 
predicts the n = 1 body mode to grow fastest at small k, as 
we have observed in Eig. (though the growth rate is only 
correct to within an 0(1) factor, as might be expected for 
such small k). 

Estimates for the surface mode frequencies can be ob¬ 
tained by equating the first and second terms in Eq. |B7| 
assuming that n is not too large. To leading order in large 
k, this balance yields 


UJ 


iqs mink 
27 4k 


(BIO) 


Thus the growth rates of the surface modes are proportional 
to qs , and their wave frequencies are typically a factor \nk/k 
smaller. 
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